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ABSTRACT 

The pressure driven, fully-developed turbulent flow of an incompressible viscous fluid in 
curved ducts of square cross-section is studied numerically by making use of a finite volume 
method. A nonlinear K - / model is used to represent the turbulence. The results for both straight 
and curved ducts are presented. For the case of fully-developed turbulent flow in straight ducts, 
the secondary flow is characterized by an eight-vortex structure for which the computed flowfield 
is shown to be in good agreement with available experimental data. The introduction of moderate 
curvature is shown to cause a substantial increase in the strength of the secondary flow and to 
change the secondary flow pattern to either a double-vortex or a four-vortex configuration. 


* This research was supported by the National Aeronautics and Space Administration under NASA Contract No. NAS1-18605 and was carried 
out while the author was in residence at ICASE, NASA Langley Research Center, Hampton, VA 23665. 
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1. Introduction 


The study of viscous flow in curved or helically coiled ducts has been of fundamental 
interest to fluid dynamicists. There are numerous applications, which include the flow through 
turbomachinery blade passages, aircraft intakes, diffusers, and heat exchangers. Some of these 
problems of practical interest involve longitudinal curvature in the geometry for which the 
associated centrifugal forces can generate a secondary flow which is normal to the main flow 
direction. Such secondary flows not only cause a reduction in the volumetric flow rate, but they 
can also cause the axial velocity field to be distorted with an outward shift of the contours of 
constant velocity. In addition, it is well known that the turbulent flow in straight noncircular 
ducts is characterized by the occurrence of secondary flows (see Gessner and Jones 1965 and 
Nakayama et al, 1983). A clear understanding of the evolution and consequences of the turbulent 
secondary flows in curved and straight ducts is, therefore, quite important from the design 
standpoint. The present study is intended to address a facet of this issue, and involves the 
computational modeling of fully developed turbulent flow in square ducts with an emphasis on 
the prediction of secondary flows. 

The results available in the literature on turbulent flow related to the work to be presented 
herein can be classified into two broad categories: fully-developed flow in straight ducts and 
developing flow in curved ducts. The case of fully-developed turbulent flow in straight ducts of 
square cross section was studied by Gessner and Jones (1965), Melling and Whitelaw (1976), 
and Nakayama, et al (1983), among others. Gessner and Jones (1965) conducted a series of 
experiments using hotwire anemometry to analyze fully-developed turbulent flow in a square 
duct at a Reynolds number of 150,000. They also carried out computations by a finite difference 
method with an algebraic stress model to predict qualitatively the major feature of the flowfield, 
namely, the eight- vortex secondary flow structure. Melling and Whitelaw (1976) performed 
detailed experiments for fully-developed flow using laser-doppler anemometry, and were the 
first to describe the axial velocity field and the Reynolds stress distribution in detail. Nakayama, 
et al (1983), on the other hand, analyzed the fully-developed flowfield in ducts of rectangular and 
trapezoidal cross-sections computationally using a finite-difference method based on the algebraic 
turbulence stress model of Launder and Ying (1972). They were able to obtain a flowfield in 
good agreement with the available experimental measurements for a number of selected cross- 
sections. Improved calculations were conducted by Gessner and Po (1977) and DeMuren and 
Rodi (1984), using the nonlinear algebraic stress model of Rodi. 
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The computational analysis of developing turbulent flow in curved square ducts has been 
conducted by Pratap and Spalding (1975), and by Humphrey, et al (1981). The work by Pratap 
and Spalding (1975) consisted of the solution of the three-dimensional time averaged Navier- 
Stokes equations incorporating a two-equation turbulence model based on a first order curvature 
ratio effect. Their study, which was conducted for a curvature ratio of 4.13 and a Reynolds 
number of 70,600, showed good agreement with the experimental results in the region near the 
entrance. However, in the fully developed region, they observed considerable discrepancies with 
the experimental results. Humphrey, et al (1981) also analyzed the developing turbulent flow in 
a square duct, but at a curvature ratio of 4.6 and a Reynolds number of 40,000. They solved the 
three-dimensional time-averaged Navier-Stokes equations with a two-equation turbulence model, 
by using a finite-difference method, and compared the results with their experimental findings 
obtained by laser-doppler anemometry. They concluded that in spite of the complex mean flow 
and Reynolds stress distributions, the cross-stream flow is controlled mainly by the centrifugal 
force and radial pressure gradient imbalance; consequently, the calculated mean velocity results 
are not strongly dependent on the turbulence model. 

The present study is motivated by the lack of a detailed analysis and prediction capability 
for turbulent secondary flows in curved ducts. Furthermore, several of the commonly used 
turbulence models (in particular, the standard K - e model) are incapable of predicting the 
development of secondary flows in straight ducts of noncircular cross-section (see Speziale 1987). 
In the computational analysis to be presented, a finite-volume algorithm suitable for handling 
cylindrical geometry is implemented to analyze fully-developed turbulent flows in straight as well 
as curved square ducts using a recently developed nonlinear K - / turbulence model (Speziale 
1987). In the following sections the governing equations, the development of the turbulence 
model, and the numerical procedure will be described in detail followed by a discussion of the 
results and conclusions. 

2. The Physical Problem and Method of Solution 

The problem to be considered consists of the turbulent flow of an incompressible viscous 
fluid in a curved duct of square cross-section. Flow in helically coiled ducts may also be analyzed 
in the same manner so long as the ratio of torsion to curvature remains small. The physical 
configuration and the coordinate system used are shown in Figure 1. 

The flow is generated by a constant azimuthal pressure gradient, G = — (see Figure 

1). It is assumed that the flow is fully-developed so that all flow variables are independent of the 
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azimuthal coordinate 0. The fully developed mean velocity field is three-dimensional, i.e. the 
mean (time-averaged) velocity vector v in the cylindrical coordinate system employed is of the 
form v = u (r,z) e r + v (r, z) e z + w (r, z) e 9 where e r , e z , ande^ denote unit vectors 
in the r,z and 6 directions, respectively. Here, u and v represent the secondary flow while w 
denotes the primary flow. The governing equations consist of the time-averaged equations for 
conservation of mass and momentum, and may be expressed in the following form 
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where G = — £ ^ is the azimuthal mean pressure gradient which is held constant, p is the fluid 
density, and v is the kinematic viscosity of the fluid which can be neglected for high Reynolds 
number turbulent flows. The components of the Reynolds stress tensor r tJ appearing in (3) can 
be obtained through various modeling techniques such as algebraic, one-equation, two-equation, 
and second-order closure models (Launder and Spalding 1972 and Lumley 1 ( >78). Recently, a 
nonlinear two-equation model of the K — l and K — e type was developed b> Speziale (1987). 
This model yields more accurate predictions for normal Reynolds stress anisotropies allowing 
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for the calculation of turbulent secondary flows. The nonlinear K — l model takes the form: 
Tij = —pK6ij -f- pK f l Sij + Co pi ^ Si m S m j ^ S mn S mn 8ij ^ 

( O 1 o 

Sij ~ Smm 8\j 


where 



are the mean rate of strain tensor and its frame-indifferent Oldroyd derivative, respectively; K 
is the turbulent kinetic energy, and / is the turbulent length scale. Cd and Ce are dimensionless 
constants that each assume the value of 1.68 which was obtained by Speziale (1987) from 
correlations with turbulent channel flow data. The turbulent length scale l can be prescribed 
algebraically by empirical means, or can be tied to the turbulent kinetic energy K and dissipation 
rate e through the relation 

r3/2 

l = 2 Cn— ( 8 ) 

where C^ is a dimensionless constant which is typically taken to be 0.09. This forms the basis 
for the nonlinear K - e model for which (5) - (7) are supplemented with modeled transport 
equations for K and e that take the form 



where cr*, a € , C t \ and C f 2 are dimensionless constants that are usually taken to be 1.0, 1.3, 1.44, 
and 1.92, respectively. The standard K - e model is obtained in the limit as Cjy, Ce — i ► 0. 

Weaknesses in the performance of the modeled transport equations (9) - (10) for K and e have 
been pointed out numerous times in the literature for problems involving swirl and streamline 
curvature (see Pope 1978 and Reynolds 1987). However, several attempts at developing 
improved versions of these modeled transport equations has not met with much success. It 
is a difficult problem that requires a substantial research effort to fully resolve. In this paper, we 
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want to examine the efficacy of the nonlinear Reynolds stress correction to the K - / and K - t 
models given by equation (5), for problems involving secondary flows with streamline curvature. 
Hence, we will specify K and / empirically, based on experimental data for turbulent flow in 
rectangular channels, in order to determine the predictive capability of the nonlinear correction to 
the Reynolds stress independent of the deficiences in the modeled transport equations for K and e. 

For fully-developed turbulent flow in a curved duct (using the cylindrical coordinate system 
shown in Figure 1), the components of the Reynolds stress tensor corresponding to the nonlinear 
K - / model (5) can be approximated as follows: 
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(see Hut 1988 for more details). In deriving (11)-(16), terms that are quadratic in the secondary 
flow velocity components u, v have been neglected since they are small (i.e., || u, v ||/|| w || <0.1 
for the computations to be presented herein where || • || denotes the maximum norm). As alluded 
to earlier, the turbulent kinetic energy K and length scale / will be specified empirically based 


5 


on the experimental data of Laufer (1951) for turbulent channel flow (i.e., for turbulent flow 
in a large-aspect-ratio rectangular duct). It will now be shown that the data of Laufer can be 
represented by the power law 
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in the interior of the duct where Uo is the centerline mean velocity, d is the half-width of the 
duct, and a *, <n, by_, and b[ are dimensionless constants. In (17) - (18), S is the mean strain 
rate defined by 
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In Figures 2-3, the turbulent kinetic energy and length scale are shown as a function of the 
mean strain rate for three different Reynolds numbers. For the following choice of empirical 
constants: 

{ 0, Sd/Uo < 0.06 

0.43, 0.06 < Sd/Uo < 0.3 (20) 

0, Sd/Uo > 0.3 


{ 0.032, Sd/Uo < 0.06 

0.11, 0.06 < Sd/Uo < 0.3 (21) 

0.064, Sd/Uo > 0.3 


{ 0, Sd/Uo < 0.04 
-0.33, 0.04 < Sd/Uo < 0.25 
-0.90, Sd/Uo >0.25 


( 22 ) 


bi 


' 0.18, Sd/Uo < 0.04 
< 0.063, 0.04 < Sd/Uo < 0.25 
k 0.028, Sd/Uo > 0.25 


(23) 


the power laws (17) - (18) do a reasonably good job in collapsing the experimental data for 
a variety of Reynolds numbers as shown in Figures 2-3. These power laws are reminiscent 
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of the ones used in viscoelastic flows; of course, the qualitative similarities between the mean 
turbulent flow of a Newtonian fluid and the laminar flow of a non-Newtonian fluid have long 
been recognized (see Lumley 1970). They have the advantage of allowing for a substantial 
reduction in the level of computation since separate transport equations for K and e do not have 
to be solved. Furthermore, they provide a more accurate measure of K and / for straight ducts 
than that which can be obtained from modeled transport equations for K and e (this allows 
for the study of the performance of the nonlinear Reynolds stress model in isolation from the 
complicating factor of defects in the K and e modeled transport equations). Although these 
power laws become less accurate for curved duct flows, the errors introduced are moderately 
small for mild curvatures and are no worse than those that arise from the standard modeled 
transport equations for K and e. 

The boundary conditions for the secondary flow velocities u and v are set equal to zero 
at the duct walls. In order to avoid resolving the very steep gradients of the azimuthal mean 
velocity w near the wall for turbulent flows, wall functions are used. These wall functions 
are based conventionally on a production equals dissipation equilibrium hypothesis and the law 
of the wall for the mean velocity profile, (c.f. Amano 1984 for more details). The boundary 
conditions for the pressure are obtained in the usual fashion (c.f. Patankar 1980). The mean 
strain rate is approximated by 



which is obtained by neglecting quadratic terms in the mean secondary flow. 

A finite volume scheme is used to solve the Reynolds equations (1) - (4), with the nonlinear 
K - / model given by equations (11) - (16). The method of solution closely follows that outlined 
by Patankar (1980) as modified for the cylindrical coordinate system which is used in conjunction 
with the curved duct geometry. In this procedure, the physical domain is discretized into a finite 
number of computational cells (see Figure 4). At the centroid of each (i.e., at point p), variables 
such as the pressure and azimuthal component of the mean velocity are defined; the components 
of the mean velocity responsible for transport in the cross-sectional planes are defined at the cell 
boundaries. The difference equations for all of the variables are then obtained by integration over 
a control volume. A detailed description of the procedure used for obtaining these equations 
may be found in Hur (1988). 
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In the present work, the system of algebraic equations which result from the differencing 
procedure used for the governing equations (1) - (4) are solved by a successive line under 
relaxation (SLUR) procedure with the repeated use of the tridiagonal matrix algorithm (Isaacson 
& Keller 1970). The details of the solution procedure and the algorithm development are given 
elsewhere (Hur 1988). The following steps constitute only a brief summary of the technique: 


(a) A set of pseudo-velocity components (for the secondary flow velocity) are obtained from the 
discretized momentum equations by assuming uniform pressure in the computational domain. 

(b) These pseudo-velocities are then used to obtain the pressure by solving the discretized 
equation for the conservation of mass with the SLUR method. 

(c) Based on the pressure obtained in (b), the secondary flow velocity components are obtained 
from the discretized momentum equations with the SLUR method. 

(d) The correction for the velocity components are obtained by evaluating the correction for the 
mass flux in each cell. 

(e) Using the corrected velocities from (d), the axial velocity is obtained by the SLUR method. 

(f) The values of kinetic energy and length scale are then computed using the corrected velocity 
field. 

(g) The procedure is repeated with the updated values of the variables until adequate convergence 
is obtained. 


In addition, a time-averaged stream function ip for the secondary flow can be defined in 
the following manner 


u 


1 dip 
r dz 


(25) 


_ 1 dip 

v = 

r or 
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and can be obtained from the mean velocity field for analyzing the secondary flowfield. In the 
following section the results obtained by the computational technique outlined herein is discussed 
and compared with available experimental and computational findings. 

3. Results and Discussion 


The model predictions for turbulent flow in a straight duct will be analyzed first, followed 
by flow in curved ducts. The experimental investigations on turbulent duct flow include the 
works of Gessner & Jones (1965), Brundrett & Baines (1964), Launder & Ying (1972), and 
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Melling & Whitelaw (1976) among others. In the present work, computations were performed 
at a Reynolds number of about 42,000 in ducts of square cross-section. This was done to 
facilitate comparison with the results of Melling & Whitelaw (1976), since their experiments 
were performed at the same Reynolds number and include laser anemometer measurements of 
the flow pattern with a detailed documentation of the Reynolds stress components. It should be 
noted here that Melling & Whitelaw (1976) did not evaluate the secondary flow stream function; 
these are obtained from other sources, although at different Reynolds numbers (e.g., Gessner & 
Jones 1965; Nakayama, et al 1983). 

The time-averaged secondary flow streamlines for the straight duct of square cross-section is 
shown in Figure 5. It should be noted that the flow is outward towards the comer and returns to 
the center along the wall thus forming eight counter-rotating vortices. Qualitative comparisons 
with the experimental results of Gessner & Jones (1965) for a Reynolds number of 150,000 
and the computational results of Nakayama, et al (1983) for a Reynolds number of 83,000 are 
possible as shown in Figure 5. The computations by Nakayama, et al (1983) are based on the 
algebraic stress model developed by Launder & Ying (1972). It can be seen here that the present 
model agrees with the experimental results of Gessner and Jones (1965) in the sense that it also 
exhibits weaker secondary flow velocities near the center of the duct. 

Due to the secondary flow associated with the turbulence, the axial velocity profile for w 
is distorted, as shown in Figure 6. Comparison with the experiments of Melling & Whitelaw 
(1976) is also shown in Figure 6 and the results indicate qualitatively good agreement. The 
computed values of the Reynolds shear stress will now be compared with the experiments. The 
results for the two shear stress components, r xz and r yz (where x and y are the horizontal and 
vertical coordinates in the plane of the secondary flow and z is the coordinate normal to this 
plane) are shown in Figures 7-8. As can be seen from these results, the computed values of 
t xz and r yz are in good agreement with the experimental findings. 

The computations were then extended to flow in curved ducts. Figures 9(a) - (d) show the 
evolution of the flow field as the curvature ratio C r (i.e., the ratio of the radius of curvature to the 
duct width) increases for a Reynolds number of about 42,000. In the case of curved ducts, the 
interaction between the centrifugal force and the force induced by the normal Reynolds stress 
differences (due to the anisotropic turbulence) characterizes the secondary flow. However, as can 
be seen from Figures 9(a) - (b), even for very loosely coiled ducts (i.e., ducts of large curvature 
ratio) the fully developed flow field is characterized by a predominately double vortex structure 
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which is representative of a flow field dominated by centrifugal effects. There is, however, a 
pair of weak counter-rotating vortices near the outer wall of the duct. When the curvature ratio 
is decreased (i.e., as the duct is coiled tighter), the centrifugal effects gain further dominance, 
and the vortices near the outer wall and, hence, the secondary flow gain strength and assume a 
four-vortex structure as shown in Figures 9(c) - (d). The effect of an increase in the secondary 
flow intensity on the azimuthal velocity profiles can also be observed in Figures 9(a) - (d). As 
can be seen from the constant velocity lines, there is a substantial outward shift of the azimuthal 
velocity field with a decrease in the curvature ratio. 

In Figure 10, the associated increase in the secondary flow intensity with a decrease in 
the curvature ratio (i.e., increase in duct curvature) is shown. From this figure, it is clear that 
curvature can increase the secondary flow intensity (defined as \u,v,\ max /w max ) from 0.02 
to 0.10. 

4. Conclusion 

A numerical study of turbulent secondary flows in straight and curved ducts of square 
cross section has been conducted using the nonlinear K - / model of Speziale (1987). The 
results are shown to compare well with the detailed laser-Doppler anemometry measurements of 
Melling & Whitelaw (1976) for fully-developed turbulent flow in straight ducts of square cross- 
section. Here, the model predicts an eight vortex secondary flow consistent with experimental 
observations. For curved ducts, the model predicts a double vortex secondary flow for mild to 
moderate curvature and a four-vortex secondary flow for moderate to strong curvature. These 
regimes are analogous to those that are observed in laminar curved duct flows (c.f. Cheng, Lin 
and Ou 1976) with one major exception: for the laminar case there would be no secondary flows 
in straight ducts (i.e., in the limit as C r — ► oo). 

In our opinion, the calculations presented in this study indicate that the nonlinear K - / 
and K - e model has the promise to yield more accurate predictions for curved turbulent flows 
than the more standardly used eddy viscosity models. This nonlinear two-equation model could 
provide a useful alternative to a second-order closure model for those applications where savings 
in computational expense is a high priority. Future research should be directed toward the 
generation of a full nonlinear K - e model solution to the curved flow problem. Invariably, this 
will require a careful examination of the modeled dissipation rate transport equation to properly 
account for curvature effects. With such improvements, the nonlinear K - e model could have a 
variety of important technological applications to turbulent flows involving streamline curvature. 
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Figure 6 
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Comparison of Reynolds Stress (T^) Contours at Re « 42000: 
Nonlinear K -/ Model; — Experiments of Melling & Whitelaw 

Figure 7 
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Comparison of Reynolds Stress (T yz ) Contours at Re * 42000: 
Nonlinear K -/ Model; - - - Experiments of Melling & Whitelaw 

Figure 8 
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1 Figure 10 
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